#include <complex>
#include "para.h"
#define C 3.0e8 //光速
#define N 12581 //所有数组的长度
#define Q 1.6e-19   //元电荷电量
#define EPSILON_0 8.854e-12 //真空介电常数
#define um 1.0e-6 //微米
#define nm 1.0e-9 //纳米
#define INT_N1 100
#define INT_N2 100
#define NX 1000
#define NZ 1000

using std::complex;

//返回单个域的相对介电常数
double d_perm(Para para, int j);
//得到所有域的相对介电常数取值数组
void d_perm(Para para, double _epsilon[], int n);



//返回单个域内的kz取值
double kz(double _epsilon[], double _omega, double _theta, int j);
//得到所有域的kz的取值数组
void kz(double _epsilon[], double _omega, double _theta, double k[], int n);



//计算单个界面的反射系数
double R(double k[], double _epsilon[], int dir, int j);
//得到每个界面的前向反射系数
void R1(double k[], double _epsilon[], double Rf[], int n);
//得到每个界面的后向反射系数
void R2(double k[], double _epsilon[], double Rb[], int n);


//计算单个界面的透射系数
double T(double k[], double _epsilon[], int dir, int j);
//得到每个界面前向透射系数
void T1(double k[], double _epsilon[], double Tf[], int n);
//得到每个界面后向反射系数
void T2(double k[], double _epsilon[], double Tb[], int n);


//计算每个界面的全局前向反射系数
void GR1(Para para, const double k[], const double Tf[], const double Rf[], const double Tb[], const double Rb[], complex<double> GRf[]);
//计算每个界面的全局后向反射系数
void GR2(Para para, const double k[], const double Tf[], const double Rf[], const double Tb[], const double Rb[], complex<double> GRb[]);


//计算因子N
void facN(Para para, const double k[], const double Tf[], const double Rf[], const complex<double> GRb[], complex<double> _N[]);
//计算因子O
void facO(Para para, const double k[], const double Tb[], const double Rb[], const complex<double> GRf[], complex<double> _O[]);



//计算每层向前辐射的振幅
void source1(const double _omega, const double par_v, Para para, const double k[], const double _epsilon[], complex<double> sf[]);
//计算每层向后辐射的振幅
void source2(const double _omega, const double par_v, Para para, const double k[], const double _epsilon[], complex<double> sb[]);


//计算前向出射的幅度
complex<double> amp_f(const int n, const complex<double> sf[], const complex<double> sb[], const complex<double> GRb[], const complex<double> _m_N[]);
//计算后向出射的幅度
complex<double> amp_b(const int n, const complex<double> sf[], const complex<double> sb[], const complex<double> GRf[], const complex<double> _m_O[]);

//计算因子I_j
complex<double> facI(Para para, int j, const complex<double> GRf[], const complex<double> GRb[], const double k[]);
//计算因子I_J数组
void facI(Para para, const int n, const complex<double> GRf[], const complex<double> GRb[], const double k[], complex<double> fI[]);

//计算“计算域j内电场强度分布需要用到的同因子N相关的连乘数组”
void relatedN(int j, const int n, const complex<double> _N[], complex<double> _rN[]);
//计算“计算域j内电场强度分布需要用到的同因子N相关的连乘数组”
void relatedO(int j, const int n, const complex<double> _O[], complex<double> _rO[]);

//计算域j内的电场强度分布
complex<double> Ezf(int j, const int n, const complex<double> sf[], const complex<double> sb[], const complex<double> GRf[], const complex<double> GRb[], const complex<double> _N[], const complex<double> _O[], const complex<double> fI[]);
complex<double> Ezb(int j, const int n, const complex<double> sf[], const complex<double> sb[], const complex<double> GRf[], const complex<double> GRb[], const complex<double> _N[], const complex<double> _O[], const complex<double> fI[]);
//计算积分
double my_int(double value[INT_N1][INT_N2], int row, int col, double dx1, double dx2);

